/*
  Copyright 2010-2011 Patrik Jonsson, sunrise@familjenjonsson.org

  This file is part of Sunrise.

  Sunrise is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 3 of the License, or
  (at your option) any later version.

  Sunrise is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Contains classes related to emission of rays from Arepo cells.

// $Id$

#ifndef __arepo_emission__
#define __arepo_emission__

#include "emission.h"
#include "arepo_grid.h"

namespace mcrx {
  template <typename, typename> class arepo_cell_emission;
}

#ifdef WITH_AREPO

/** Represents isotropic emission from an Arepo Voronoi cell. The
    position within the Voronoi cell is drawn randomly from a
    distribution proportional to the density within the cell using
    rejection sampling. The set up of the bounding box information is
    a bit problematic as this requires the circumsphere centers of the
    Delaunay tessellation, which is not stored indexed by the
    cells. For this reason, after creating all the cell_emission
    objects, a final pass is done which loops over the Delaunay
    tetrahedrons and updates the bounding box information for the
    cells associated with that tetrahedron. This is done by the static
    post_setup() function, which must be called after the grid has
    been constructed. This is unfortunate, but the alternative is to
    loop over the DT array once per created cell, which would be
    prohibitively expensive.  */
template <typename chromatic_policy, 
	  typename rng_policy_type>
class mcrx::arepo_cell_emission :
  public isotropic_emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename isotropic_emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename isotropic_emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef rng_policy_type T_rng_policy;

  /** To avoid circularly dependent template parameters, it is
  hardcoded that the grid cells contain this class.  This precludes
  making grids which contain objects which in turn contains
  grid_cell_emission objects, but was chosen for simplicity and the
  lack of the current need for anything more advanced. */
  typedef arepo_grid_cell<arepo_cell_emission> T_cell;

private:
  /** True if the bounding box setup has been done. This is just used
      to check that the post_setup step has been done. */
  bool setup_done_; 

  const T_cell* cell_; ///< Pointer to the cell doing the emission.

  T_float rho_min_; ///< The mininimum density in the cell
  T_float rho_max_; ///< The maximum density in the cell

  vec3d min_; ///< The min coordinates of the cell bounding box
  vec3d max_; ///< The max coordinates of the cell bounding box

  T_lambda emission_; ///< The "thing" being emitted.

public:
  arepo_cell_emission (const T_lambda& em, const T_cell& c):
    isotropic_emission<chromatic_policy, rng_policy_type>(em), 
    setup_done_(false), cell_(&c), emission_(em), 
    rho_min_(SphP[cell_->index()].Density), rho_max_(rho_min_),
    min_(cell_->position()), max_(cell_->position())
  {
    // Note that the bounding box can not be set correctly at this
    // point. These must be updated by looping over the tetrahedrons
    // that are connected to this mesh point.
    ASSERT_ALL(em==em);
    ASSERT_ALL(em<1e300);
  };

  /** Update the bounding box information with the fact that p is a
      point on the cell face. */
  void update_bounding_box(const vec3d& p) {
    min_=where(p<min_, p, min_); max_=where(p>max_, p, max_);
    // update density min and max values with the values at this location
    const T_float rho = arepo::cell_density(cell_->index(), p);
    assert(rho>=0);
    rho_min_=std::min(rho_min_, rho);
    rho_max_=std::max(rho_max_, rho);
 };

  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy&); 
  T_lambda zero_lambda() const {return T_lambda(0.*emission_);}; 
  std::auto_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return std::auto_ptr<emission<chromatic_policy, rng_policy_type> >(new arepo_cell_emission<chromatic_policy, rng_policy_type>(*this));};

  // copy constructor and assignment operator are implicitly generated

  const T_lambda& get_emission() const { return emission_;};

  /** This function sets the bounding boxes of the cells by looping
      over all the Arepo tetrahedrons and updating the bounding boxes
      with the positions of the circumsphere centers of tetrahedrons
      connected to the mesh point. */
  static void update_bounding_boxes(const arepo_grid<arepo_cell_emission>& g);

  /** This routine uses the arepo_grid object to setup of the bounding
      box information. It MUST be called before emitting rays. */
  static void post_setup (const arepo_grid<arepo_cell_emission>& g) {
    update_bounding_boxes(g); };
};


/** Emit a ray from a location within the cell drawn based on the
    density distribution. The direction is drawn isotropically. */
template <typename chromatic_policy, 
	  typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::arepo_cell_emission<chromatic_policy, rng_policy_type>::T_lambda>
mcrx::arepo_cell_emission<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng)
{
  using namespace arepo;
  assert(setup_done_);

  // we draw a point consistent with the density distribution in the
  // cell using rejection sampling
  const int dp_idx=cell_->primary_cell();
  const int sphp_idx = sph_index(dp_idx);

  const vec3d cell_p (get_point(dp_idx));
  vec3d pos;
  int n=0;

  DEBUG(1,std::cout << "arepo_cell_emission: Drawing emission point for cell " << sphp_idx << endl;);

  while(true) {
    ++n;
    DEBUG(1,std::cout << "\tAttempt " << n << "\n";);
    assert(n<10000);

    // first we draw a point from within the bounding box
    blitz::TinyVector<T_float, 3> r;
    rng.rnd(r);
  
    // randomly draw a position
    pos = (max_-min_)*r + min_;

    // test if point is within cell
    bool inside=true;
    const std::pair<int, int> edge = arepo::midpoint_idx[dp_idx];

    for(int i=edge.second-1; i>=0; --i) {
      // c is the vector from entry point p0 to the point on the face
      const vec3d c(arepo::midpoints[edge.first+i] - pos);
    
      /* q is (half) the edge vector to the neighboring cell, which is a
	 normal vector of the plane */
      const vec3d q(arepo::midpoints[edge.first+i] - cell_p);

      // The point is on the inside of this face if c.q>0.
      const T_float cdotq = dot(c,q);
      if (cdotq<=0) {
	// point is on outside of this face, we can stop testing now
	inside=false;
	break;
      }
    }
    // if we come out here and inside is true, the point is inside.
    if(!inside) {
      // If it's outside we try a new position
      DEBUG(1,std::cout << "\tPoint " << pos << " outside cell, rejecting\n";);
      continue;
    }

    DEBUG(1,std::cout << "\tPoint " << pos << " inside cell, proceeding...\n";);

    // if we are inside, we now accept with a probability of rho(pos)/rho_max.
    const T_float r2(rng.rnd());
    if (r2*rho_max_ < rho_min_) {
      // we can accept without testing, because we know the density in
      // this spot is larger than rho_min_.
      DEBUG(1,std::cout << "\taccepting density without testing\n";);
      break;
    }

    // we could not accept without testing, so we need to calculate
    // the actual density at this point
    const T_float rho = arepo::cell_density(dp_idx, pos);
    if (r2*rho_max_ < rho) {
      // we accept this point
      DEBUG(1,std::cout << "\tDensity fraction " << r2 << " accepted.\n";);
      break;
    }
    DEBUG(1,std::cout << "\tDrawn density fraction " << r2 << " higher than " << rho/rho_max_ << ", rejecting\n";);
  }

  // we come out here with a point correctly drawn according to the
  // density distribution in the cell. Now draw an angle and we're
  // done.
  blitz::TinyVector<T_float, 2> r;
  rng.rnd(r);

  // and the emission direction (from an isotropic distribution)
  const T_float dphi = r [0]*2*constants::pi; 
  const T_float dcos_theta = 2*r [1] - 1;
  const T_float dsin_theta = sqrt (1- dcos_theta*dcos_theta);

  //and create the ray
  return ray_distribution_pair<T_lambda> 
    (boost::shared_ptr<T_ray> 
     (new T_ray (pos,
		 vec3d (dsin_theta*cos (dphi), dsin_theta*sin (dphi), 
			dcos_theta),
		 initialize(emission_, 1.))), 
     emission_,
     this->ad_);
}

template <typename chromatic_policy, 
	  typename rng_policy_type>
void
mcrx::arepo_cell_emission<chromatic_policy, rng_policy_type>::
update_bounding_boxes(const arepo_grid<arepo_cell_emission>& g)
{
  using namespace arepo;
  tetra* DT = T->DT;
  point* DP = T->DP;
  tetra_center* DTC = T->DTC;
  std::cout << "Calculating bounding boxes of cells." << std::endl;

  // Start by looping over all tetras
  for(int t=0; t<T->Ndt; ++t) {
    if( DT[t].t[0]<0 || 
	DT[t].p[0] == DPinfinity || DT[t].p[1] == DPinfinity || 
	DT[t].p[2] == DPinfinity || DT[t].p[3] == DPinfinity )
      // tetras with t[0]<0 or that contain an infinity point do not
      // have computed circumspheres. (It's fine, <0 tetras are
      // deleted and the infinity points are well outside our box.)
      continue;

    const vec3d dtc (vec3d(DTC[t].cx, DTC[t].cy, DTC[t].cz)*arepo::lcon);

    // For each tetrahedron, we loop over the four vertices
    for(int p=0; p<4; ++p) {
      // find the cell connected to this vertex
      const int dp_idx =DT[t].p[p];
      if(arepo::primary_cell(dp_idx)) {
	// primary cell, get the emission object
	const int sphp_index = sph_index(dp_idx);
	arepo_cell_emission* c = g.cell(sphp_index).data();
	// and update bounding box
	c->update_bounding_box(dtc);
	c->setup_done_=true;
      }
    }
  }
}

#endif
#endif
